package scu.maqiang.possion;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PossionFPI {

	public static void main(String[] args) {
		Mesh2T3 mesh = new Mesh2T3().square2D(100, 100, T3Type.Cross1);
		FES2T31 fs = new FES2T31(mesh);
		
		int dof = fs.GetNdof();
		SRMatrix A = new SRMatrix(dof);
		double[] coef = new double[dof];
		double[] RHS = new double[dof];
		double[] x = new double[dof];
		double[] x0 = new double[dof];
		double[][] myCoef = new double[1][];
		fs.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
		NewIterSSolver solver = new NewIterSSolver(A);
		double err = 1.0;
		while(err > 1.0e-10) {
			Arrays.setAll(coef, i -> 2.0 + x[i] + x[i] * x[i] * x[i]);
			//Arrays.setAll(coef, i -> 1.0);
			myCoef[0] = coef;
			A.clear();
			fs.assembleStiff(myCoef, BVPType.COMMON, A);
			fs.applyBCZero(A, RHS, Direct.All, BCType.MBN, 1, 2, 3, 4);
			solver.setMatrix(A);
			solver.PCGSSOR(RHS, x, 1.5, 1);
			err = MVO.L2Norm(MVO.add(x, -1, x0));
			Arrays.setAll(x0, i -> x[i]);
		}
		mesh.toTecplot("FPIResult.dat", x);
	}

}
